function a=DLSQR(x,y)
A = preparematrix(x);
[Q,R] = qr(A);
R1 = R(1:3,:);
c = Q'*y';
c = c(1:3);
a = R1 \ c;